home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / fraction < prev    next >
Text File  |  1996-04-09  |  2KB  |  68 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    convert a real number to a nearest proper fraction.
  4.  
  5. // Syntax:      fraction ( num , maxerr )
  6.  
  7. // Description:
  8.  
  9. // Fractional approximation of a real number num.  Return a row vector
  10. // [a, b, c] where a is the integral part of num, b is the numerator,
  11. // and c is the denominator.  The maximum error can be specified in
  12. // maxerr such that abs(num - a - b/c) < maxerr.  If maxerr is not
  13. // specified, then 1/64 wil be used as maxerr.
  14. //
  15. // Examples:
  16. //
  17. // > fraction(0.33)
  18. //         0          1          3
  19. // > fraction(0.33, 0.002)
  20. //         0         21         64
  21. // > fraction(0.33, 0.001)
  22. //         0         26         79
  23. // > fraction(1.33)
  24. //         1          1          3
  25. // > fraction(-1.25)
  26. //        -1         -1          4
  27. //
  28. // Reference: L. J. Plebani, "Common-fraction approximation of real
  29. //               numbers," Dr. Dobb's J., October, 1995.
  30. //
  31. // Tzong-Shuoh Yang  9/8/95   (yang@isec.com)
  32. // (I love SI much much more than English unit system)
  33. //-------------------------------------------------------------------//
  34. fraction = function(num, maxerr)
  35. {
  36.    local(num, maxerr)
  37.  
  38.    if(!exist(maxerr))
  39.    {
  40.        maxerr = 1/64;
  41.    }
  42.    num = real(num);
  43.    maxerr = abs(real(maxerr));
  44.  
  45.    if (num > 0) { s = 1; else s = -1; num = abs(num);}
  46.    i = int(num);
  47.    num = num - i;
  48.    a = 0; b = 1; c = 1; d = 1;
  49.    while(1)
  50.    {
  51.        h = a + c;
  52.        k = b + d;
  53.        e = h/k - num;
  54.        if (abs(e) < maxerr)
  55.        { 
  56.            break;
  57.        else if (e < 0)
  58.        { 
  59.            a = h;
  60.            b = k;
  61.        else
  62.            c = h;
  63.            d = k;
  64.        }}
  65.    }
  66.    return [s*i,s*h,k];
  67. };
  68.